####install/load relevant packages####

install.packages("metafor")
install.packages("ggplot2")
install.packages("rotl")
install.packages("ape")
install.packages ("lme4")

library(metafor)
library(ggplot2)
library(rotl)
library(ape)
library(lme4)

#to install the latest version of orchaRd, use the following code (source: https://github.com/daniel1noble/orchaRd)

install.packages("pacman")
pacman::p_load(devtools, tidyverse, metafor, patchwork, R.rsp, emmeans)
devtools::install_github("daniel1noble/orchaRd", force = TRUE)
library(orchaRd)

####R2 Function####

# Function for R2 provided by Shinichi
# make a version which provides both R types of R^2

R2 <- function(model){
  warning("Make sure you have the observation (effect size) level random effect as the last in the formula\n")
  
  # fixed effect variance
  fix <- var(as.numeric(as.vector(model$b) %*% t(as.matrix(model$X))))
  
  # marginal
  R2m <- fix / (fix + sum(model$sigma2))
  R2
  #Rm <- round(100*R2m, 3)
  
  # conditional
  R2c <- (fix + sum(model$sigma2) - model$sigma2[length(model$sigma2)]) / 
    (fix + sum(model$sigma2))
  
  R2s <- c(R2_marginal = R2m, R2_conditional = R2c)
  
  return(R2s)
}

####phylogenetic tree and matrix####

fishinb<-read.csv(file.choose())

fishinb$GENUS_SPECIES<-as.factor(fishinb$GENUS_SPECIES) #set as a categorical factor
fishinb$GENUS_SPECIES_corr<-as.factor(fishinb$GENUS_SPECIES_corr) #set as a categorical factor

#code below checks that all of our species can be found in the Open Tree of Life database#

species_levels <- levels(fishinb$GENUS_SPECIES)
matchnames <- tnrs_match_names(species_levels)
check <- matchnames[matchnames$number_matches>1,]
check
ott_list <- as.list(check$ott_id)
options <- lapply(ott_list, function(x) inspect(matchnames, ott_id = x))
options
choose <- c(1,
            2)
take <- list()
for(i in 1:length(choose)){
  take[[i]] <- options[[i]][choose[i],]
}
take <- take %>% bind_rows()
rawtree <- tol_induced_subtree(ott_ids = matchnames$ott_id, label_format = "name")

#code below builds the phylogenetic tree and assigns branch lengths using the Grafen method
#branch lengths initially set to 1, then made ultrametric

tips <- gsub("_", " ", rawtree$tip.label)
tips <- gsub(" *\\(.*?\\) *", "", tips)
extra <- tips[(tips %in% species_levels) == F]
miss <- species_levels[(species_levels %in% tips) == F]
tips[(tips %in% species_levels) == F] <- species_levels[(species_levels %in% tips) == F]
rawtree$tip.label <- tips
rawtree_full <- rawtree
set.seed(49)
bintree_full <- multi2di(rawtree_full, random = T) 
is.binary(bintree_full) # TRUE
phylo_branch_full <- compute.brlen(bintree_full, method = "Grafen", power = 1)
plot(phylo_branch_full,cex=0.6)
is.ultrametric(phylo_branch_full) # TRUE

#code below generates the variance-covariance matrix of phylogenetic relatedness ("phylo_cor_full") that we then use as a random effect in our multi-level meta-analytic models

phylo_cor_full <- vcv(phylo_branch_full, cor = T)
phylo_branch_full$node.label <- NULL
phylo_MCMC <- MCMCglmm::inverseA(phylo_branch_full, nodes = "ALL", scale = TRUE)$Ainv

####multi-level meta-analytic models####

#the code below sets these variables as categorical factors#

fishinb$ ..Paper_ID<-as.factor(fishinb$ ..Paper_ID)
fishinb$Dvar_ID<-as.factor(fishinb$Dvar_ID)
fishinb$Care.provided.<-as.factor(fishinb$Care.provided.)
fishinb$Sex.providing.care<-as.factor(fishinb$Sex.providing.care)
fishinb$Timing.of.care<-as.factor(fishinb$Timing.of.care)
fishinb$Fitness.parameter<-as.factor(fishinb$Fitness.parameter)
fishinb$Study.setting<-as.factor(fishinb$Study.setting)
fishinb$Life.stage<-as.factor(fishinb$Life.stage)
fishinb$Focal.sex<-as.factor(fishinb$Focal.sex)
fishinb$Stress<-as.factor(fishinb$Stress)

summary(fishinb) #check that everything looks correct

#REM1: main effects model with no moderator (testing whether there is evidence for significant inbreeding depression overall in the dataset#

REM1=rma.mv(yi=Hedges.g,V=Hedges.g.var, random = list(~1| ..Paper_ID, ~1|Shared.control, ~1|GENUS_SPECIES, ~1|GENUS_SPECIES_corr, ~1|Dvar_ID), R = list(GENUS_SPECIES_corr = phylo_cor_full), method="REML",data=fishinb)

summary(REM1) #this produces the output for the main effects model

#REM2: moderator model with a continuous variable (F-value)

REM2=rma.mv(yi=Hedges.g,V=Hedges.g.var, mod=~F.value, random = list(~1| ..Paper_ID,  ~1|Shared.control, ~1|GENUS_SPECIES, ~1|GENUS_SPECIES_corr, ~1|Dvar_ID,), R = list(GENUS_SPECIES_corr = phylo_cor_full), method="REML",data=fishinb)

summary(REM2) #this produces the output for moderator model REM2
R2(REM2) #this calculates the marginal R2, which tells us how much heterogeneity in effect sizes is explained by the moderator

#REM3a: moderator model with a categorical variable (whether or not care was provided)

REM3a=rma.mv(yi=Hedges.g,V=Hedges.g.var, mod=~Care.provided., random = list(~1| ..Paper_ID, ~1|Shared.control, ~1|GENUS_SPECIES, ~1|GENUS_SPECIES_corr, ~1|Dvar_ID), R = list(GENUS_SPECIES_corr = phylo_cor_full), method="REML",data=fishinb)

summary(REM3a)
R2(REM3a)

#REM3a tells us whether there is a significant effect of the moderator - we then need to run a no-intercept model (REM3b) with this moderator to obtain mean effect size estimates and confidence intervals for each moderator level (we do this by adding "-1" after the moderator)

REM3b=rma.mv(yi=Hedges.g,V=Hedges.g.var, mod=~(Care.provided.)-1, random = list(~1| ..Paper_ID, ~1|Dvar_ID, ~1|Shared.control, ~1|GENUS_SPECIES, ~1|GENUS_SPECIES_corr), R = list(GENUS_SPECIES_corr = phylo_cor_full), method="REML",data=fishinb)

summary(REM3b)

REM4a=rma.mv(yi=Hedges.g,V=Hedges.g.var, mod=~(Sex.providing.care), random = list(~1| ..Paper_ID, ~1|Shared.control, ~1|GENUS_SPECIES, ~1|GENUS_SPECIES_corr, ~1|Dvar_ID), R = list(GENUS_SPECIES_corr = phylo_cor_full), method="REML",data=fishinb)

summary(REM4a)
R2(REM4a)

REM4b=rma.mv(yi=Hedges.g,V=Hedges.g.var, mod=~(Sex.providing.care)-1, random = list(~1| ..Paper_ID, ~1|Dvar_ID, ~1|Shared.control, ~1|GENUS_SPECIES, ~1|GENUS_SPECIES_corr), R = list(GENUS_SPECIES_corr = phylo_cor_full), method="REML",data=fishinb)

summary(REM4b)

#the linear contrasts below are performing pairwise comparisons for the moderator levels of parental care (maternal, biparental, or no care) to determine which ones are significantly different from each other

anova(REM4b, L=c(0,1,-1)) #Sex.providing.careFemale - Sex.providing.careNone = 0 
anova(REM4b, L=c(1,-1,0)) #Sex.providing.careBoth - Sex.providing.careFemale = 0 
anova(REM4b, L=c(1,0,-1)) #Sex.providing.careBoth - Sex.providing.careNone = 0 

REM5a=rma.mv(yi=Hedges.g,V=Hedges.g.var, mod=~(Timing.of.care), random = list(~1| ..Paper_ID, ~1|Shared.control, ~1|GENUS_SPECIES, ~1|GENUS_SPECIES_corr, ~1|Dvar_ID), R = list(GENUS_SPECIES_corr = phylo_cor_full),  method="REML",data=fishinb)

summary(REM5a)
R2(REM5a)

REM5b=rma.mv(yi=Hedges.g,V=Hedges.g.var,mod=~(Timing.of.care)-1, random = list(~1| ..Paper_ID, ~1|Dvar_ID, ~1|Shared.control, ~1|GENUS_SPECIES, ~1|GENUS_SPECIES_corr), R = list(GENUS_SPECIES_corr = phylo_cor_full), method="REML",data=fishinb)

summary(REM5b)

REM6a=rma.mv(yi=Hedges.g,V=Hedges.g.var,mod=~(Fitness.parameter), random = list(~1| ..Paper_ID, ~1|Dvar_ID, ~1|Shared.control, ~1|GENUS_SPECIES, ~1|GENUS_SPECIES_corr), R = list(GENUS_SPECIES_corr = phylo_cor_full), method="REML",data=fishinb)

summary(REM6a)
R2(REM6a)

REM6b=rma.mv(yi=Hedges.g,V=Hedges.g.var,mod=~(Fitness.parameter)-1, random = list(~1| ..Paper_ID, ~1|Dvar_ID, ~1|Shared.control, ~1|GENUS_SPECIES, ~1|GENUS_SPECIES_corr), R = list(GENUS_SPECIES_corr = phylo_cor_full), method="REML",data=fishinb)

summary(REM6b)

REM7a=rma.mv(yi=Hedges.g,V=Hedges.g.var,mod=~(Study.setting), random = list(~1| ..Paper_ID, ~1|Shared.control, ~1|GENUS_SPECIES, ~1|GENUS_SPECIES_corr, ~1|Dvar_ID), R = list(GENUS_SPECIES_corr = phylo_cor_full), method="REML",data=fishinb)

summary(REM7a)
R2(REM7a)

REM7b=rma.mv(yi=Hedges.g,V=Hedges.g.var, mod=~(Study.setting)-1, random = list(~1| ..Paper_ID, ~1|Dvar_ID, ~1|Shared.control, ~1|GENUS_SPECIES, ~1|GENUS_SPECIES_corr), R = list(GENUS_SPECIES_corr = phylo_cor_full), method="REML",data=fishinb)

summary(REM7b)

REM8a=rma.mv(yi=Hedges.g,V=Hedges.g.var,mod=~(Life.stage), random = list(~1| ..Paper_ID, ~1|Dvar_ID, ~1|Shared.control, ~1|GENUS_SPECIES, ~1|GENUS_SPECIES_corr), R = list(GENUS_SPECIES_corr = phylo_cor_full), method="REML",data=fishinb)

summary(REM8a)
R2(REM8a)

REM8b=rma.mv(yi=Hedges.g,V=Hedges.g.var,mod=~(Life.stage)-1, random = list(~1| ..Paper_ID, ~1|Dvar_ID, ~1|Shared.control, ~1|GENUS_SPECIES, ~1|GENUS_SPECIES_corr), R = list(GENUS_SPECIES_corr = phylo_cor_full), method="REML",data=fishinb)

summary(REM8b)

#REM9a and REM9b filters out effect sizes where the focal sex is categorised as both, so that we can directly compare the magnitude of inbreeding depression in males vs females

REM9a=rma.mv(yi=Hedges.g,V=Hedges.g.var, mod=~(Focal.sex), random = list(~1| ..Paper_ID, ~1|Dvar_ID, ~1|Shared.control, ~1|GENUS_SPECIES, ~1|GENUS_SPECIES_corr), R = list(GENUS_SPECIES_corr = phylo_cor_full), method="REML",
             data=fishinb %>% filter(fishinb$Focal.sex != "Both"))

summary(REM9a)
R2(REM9a)

REM9b=rma.mv(yi=Hedges.g,V=Hedges.g.var,mod=~(Focal.sex)-1, random = list(~1| ..Paper_ID, ~1|Dvar_ID, ~1|Shared.control, ~1|GENUS_SPECIES, ~1|GENUS_SPECIES_corr), R = list(GENUS_SPECIES_corr = phylo_cor_full), method="REML",
             data=fishinb %>% filter(fishinb$Focal.sex != "Both"))

summary(REM9b)

REM10a=rma.mv(yi=Hedges.g,V=Hedges.g.var,mod=~(Stress), random = list(~1| ..Paper_ID, ~1|Dvar_ID, ~1|Shared.control, ~1|GENUS_SPECIES, ~1|GENUS_SPECIES_corr), R = list(GENUS_SPECIES_corr = phylo_cor_full), method="REML",data=fishinb)

summary(REM10a)
R2(REM10a)

REM10b=rma.mv(yi=Hedges.g,V=Hedges.g.var,mod=~(Stress)-1, random = list(~1| ..Paper_ID, ~1|Dvar_ID, ~1|Shared.control, ~1|GENUS_SPECIES, ~1|GENUS_SPECIES_corr), R = list(GENUS_SPECIES_corr = phylo_cor_full), method="REML",data=fishinb)

summary(REM10b)

####publication bias####

REMpubbias=rma.mv(yi=Hedges.g,V=Hedges.g.var,mod=~Year, random = list(~1| ..Paper_ID, ~1|Dvar_ID, ~1|Shared.control, ~1|GENUS_SPECIES, ~1|GENUS_SPECIES_corr), R = list(GENUS_SPECIES_corr = phylo_cor_full), 
                  method="REML",data=fishinb)

summary(REMpubbias)

REM_pubbias=rma(yi=Hedges.g, vi=Hedges.g.var,method="REML",data=fishinb)
summary(REM_pubbias)

pubbiasplot <- funnel(REM_pubbias,cexr.main=1.25,cex.lab=1,cex.axis=1,cex=1)
regtest(REM_pubbias) #Test for funnel plot asymmetry

####POWER ANALYSIS####

devtools::install_github("dsquintana/metameta")
install.packages("metameta")
library(metameta)

install.packages("remotes")
remotes::install_github("dsquintana/metameta")

dat_pp<-read.csv(file.choose())

dat_pp$yi <- dat_pp$Hedges.g
dat_pp$sei <- dat_pp$Hedges.g.var

power_pp <- mapower_se(dat = dat_pp, 
                        observed_es = 0.20, 
                        name = "P & P")
power_pp_dat <- power_pp$dat
power_pp_dat

power_pp <- mapower_se(dat = dat_pp, observed_es = 0.20, name = "P & P")
power_med_pp <- power_pp$power_median_dat

list_power <- list(power_med_pp)
